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Abstract 

We study quantum phase transition from the superfluid to a Mott insulator in optical lattices us- 
■ ing a Bose-Hubbard Hamiltonian. For this purpose we have developed a field theoretical approach 
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^2 ! density as well as to the superfiuid fraction in cubic optical lattices. Using present approach the 



condensate fraction and ground state energy are calculated as functions of the s-wave scattering 
length. In contrast to the Bogoliubov model, which is technically speaking a one-loop approxi- 
mation, we carry the calculation up to two loops, and improve the result further by variational 
perturbation theory. The result suggests that the quantum phase transition exists. 
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I. INTRODUCTION 



Optical lattices are known as the gases of ultracold atoms trapped in periodic potentials 
created by standing waves of laser light. The actuality of experimental and theoretical 
investigations of these artificial crystals bound by light can be justified by following two 
factors [l|: 

1) Neutral atoms in these optical lattices have a number of affective futures that make 
them interesting candidates for the realization of a quantum computer [2]. 

2) They may be used to stimulate various lattice models of fundamental importance to 
condensed matter physics to study in a controlled way in solid-state physics, since one is 
able to finely tune the properties and geometry of the lattices. In particular, it is possible 
to control the Hamiltonian parameters and study various regimes of interest. Similarly to 
the ordinary Bose - Einstein Condensation (BEC) of gases, the quantum phase transitions 
in optical lattices were first predicted theoretically |3| and have recently been observed 
experimentally 14 1. 

Most of the theoretical investigations are based on Bose-Hubbard Hamiltonian: 

Tj N, Ns 

h = -jy: 4ci + ^T.^4ci£i + T.i^i-Mci (1) 

<ij> i i 

where Ci^ and q are the bosonic creation and annihilation operators on the site i; the sum 
over < i, j > includes only pairs of nearest neighbors; J is the hopping amplitude, which is 
responsible for the tunneling of an atom from one site to another neighboring site; U is the 
on site repulsion energy; Ng - number of sites. Presently it is well established that at very 
low temperature (T — > 0) a system of bosons described by the Hamiltonian ([1]) could be on 
superfluid (SF) or in Mott insulator (MI) phase. Clearly there would be a quantum phase 
transition between these two phase depending on parameters U and J. Particularly, when 
the hopping term is dominated, U/J'^ 1, the system prefers to be in the SF phase. On the 
other hand when the repulsion prevails the kinetic term, U/J^ 1, the system would be in 
MI phase where each atoms is absolutely localized near a site. 

Clearly the superfluid phase may consist not only of condensed particles with a number 
A'^o, but also of A^^i uncondensed ones, whose sum A^^q + A^^i = is the total number of 
particles. The critical interaction strength Kent = (f^/-^)crit ~ 29.34 and Kent = 3.6, for 
D = 3 and D = 1 respectively, of the quantum phase SF MI transition estimated 
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by Monte Carlo calculations at filling factor z/ = 1 is in good agreement with the 

experimental data. 

To make easier further reading we clarify some specific features of these two phases: 
SF phases is characterized by long-range correlation, a continuous (gapless) excitation spec- 
trum and a finite compressibility. Since there exists a condensate with a number of particles 
A^o 7^ 0, the gauge symmetry is spontaneously broken in accordance with Bogoliubov and 
Ginibre theorems. In contrast, in the Mott insulator phase, there is no long-range correlation 
neither breaking of gauge symmetry. The excitation spectrum has a gap and the system is 
incompressible, since there is a fixed number of atoms per-site. This new state of matter 
can survive only at zero temperature and integer filling factor u. 

It is interesting to note that there are two kinds of experiments observing above quantum 



phase transition, depending on the starting point. In the experiments by Greiner et al [4] 
one first creates a BEG in a conventional harmonic trap and then adiabatically adds the 
periodic optical potential. In the second method, pioneered by the Florence group [3] one 
uses a conventional protocol for evaporative cooling in a magnetic trap down to temperatures 
just above the threshold for BEG. At this point the optical lattice potential is switched on 
and evaporative cooling continues. In this way, the system condenses directly into a ground 
state of the harmonic plus periodic potential. It seems to be that the first method is good 
to observe SF— )-MI while the second one is good for MI— t-SF transitions. 

Similarly, most of theoretical approaches can be divided into two classes: SF— t-MI and 
MI— j-SF ones. The latter are based on the Ginzburg - Landau theory as describes for 
instance in Ref. jsl. They are well suited to analyze the time-of light pictures and the 
resulting visibility at zero and finite temperatures. In the former class (SF— J-MI) one uses a 



perturbative scheme jol within a decoupling (or single site) approximation due to Gutzwiller. 



This variational appoach which was first proposed for a fermion system 10|, and further 



developed for bosons in Refs. UJ, ll2| , has the following drawbacks 
of Sec. IV): 



13-1 (see also last hues 



The mean field Hamiltonian which features single boson terms does not conserve the 

Tunneling of uncondensed atoms is neglected; 

The critical value Kent does not depend on the lattice dimension. 



Nevertheless, the prediction of decouphng approximation for Kent = 34.98 at filhng factor 
= 1 is in agreement with the well established value given above. Some years ago an 
application of the Hartree-Fock- Popov approximation (which is widely used to study BEC 



of atomic gases and even triplons 15|) to optical lattices was presented by Stoof et al. 16 | 
Studying the dependence of the condensate number Nq on k = U/J, i.e. Nq{U / J) they 
observed that A^o never reaches zero for finite values of k, implying that this approximation 
is unable to predict a possible phase transition to a Mott-insulator phase. Moreover, a 
Hartree-Fock-Bogoliubov (HFB) approximation applied to the Bose-Hubbard Hamiltonian 
gives no quantum phase transition for optical lattices Hence we find it interesting to 
study the possibility of such a transition if we go beyond these approximations. 

In the present work we shall investigate BEC in optical lattices by applying a two-loop 
approximation and treating the result by variational perturbation theory (VPT) [18]. It 
will be shown that, while the ground state energy is rather sensitive to the filling factor in 
commensurate situations, this is not so for arbitrary condensate fractions hq = No/N. We 
find that hq goes to zero at k ~ 6 ^ 6.5 for z/ = 1, 2, 3 in D = 3 dimensions. In D = 1 
dimension, this happens at k ~ 4. 

The plan of this paper is as follows. In Sec. 11 the basic equations in functional formalism 
for Bose-Hubbard Hamiltonian are formulated. In Sec. HI we derive explicit expressions for 
the effective potential in two-loop order. In Sec. IV we obtain condensate fraction vs input 
parameters f/, J, z/ . The quantum corrections to the energy of the system is discussed in Sec. 
V. In Sec. VI we present numerical results and discussions. The last Sec. VII summaries 
our results. 

II. THE ACTION AND PROPAGATORS IN BOSE-HUBBARD MODEL 

The action at zero temperature, (T = 0) that describes a gas of atoms in a periodic 
potential is given by 

/■ r 

S{(p\(p) = J dtdic ipUdtip + (p^—(p + iJ,(p^(p-Vc^ti^)(p^(p 

-^J y?^(x)(/?t(x')V(x - x')</?(x)^(x')rftc/xrfx' (2) 
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where the isotropic optical lattice potential is described by j4| 

^ext(x) = yoEsin2(^) (3) 
with A the wave length of the laser light. The lattice points lie at the positions [l^ 

Xi = i a, (4) 

where a is the lattice spacing, and 

i = {ii,i2, ■■■,id) (5) 

are integer-valued vectors. It can be shown j^, [isl that the Wannier representation of the 
Hamiltonian corresponding to the action ([2]) is equivalent to well known Bose-Hubbard 
model (d]). 

The on-site energy, ^i, the amplitude of hopping - J and on-site interaction strength U 
are related to V"ext(x) and V{'x — x') as follows: 



2m 



+ Kxt(x) ^u;o(x-Xi) (6) 



£i = y" d^uliy. - Xi) |- 

/I Ti^W^ I 
(ixa;(5(x - Xi) I — ^— + Kxtlx) \ c<;o(x - xj) (7) 

U = J J c?x'(X'd(x — Xi)u;d(x — Xi)l^(x — x')(X'o(x' — Xi)(jJo(x' — Xi) (8) 

where a;„(x) are Wannier functions. In the tight-binding limit and pseudopotential approx- 
imation, \^(x — x') = 47ra5(x — x')/m the equations ([7]), ([8]) are simplified as: 



u = ^ (10) 



where Er = 27r^/mA^, a is the s-wave scattering length, and / = yl/mu = (i?r/Vo)^^^A/47r 
is the harmonic oscillator length. 

In terms of parameters J and U the action can be rewritten as follows: 

S{v\^) = / dtl ^(p\xi,t)[idt + fx\(fi{xi,t) + J (/?^(xi,t)(/?(xj,t) 

I- i <i,j> 

-^Yl V'H^i' ^)v\^u i)v5(xi, t)v?(xi, t) \ (11) 
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The Krand-canonical partition function Z, and the effective potential at zero temperature, 
V, can be found as [20|: 



Z 
V 



InZ 



(12) 
(13) 



where J dt = T is the total time interval. Note that, in accordance with the background 
field method 2l|, which will be used below, in evaluation of the effective potential only 



connected single - particle irreducible Feynman diagrams should be included. The ground 
state expectation value of an operator A{ip'^ , ip) can be expressed as a functional integral: 



(14) 



At zero temperature the system could undergo into BEC state. The necessary and sufficient 
condition for Bose-Einstein condensation is the spontaneous gauge-symmetry breaking which 
is established by Bogoliubov shift [l3|: 



V9(xi, t) = y/un^ + if{^i,t) 



(15) 



where i' = N/Ng- filling factor, and the condensate fraction, no = Nq/N, is constant for 
regular lattice without magnetic trap. 

Substituting (fT5|) into (fTTl) and parameterizing quantum field <^(xi,t) in terms of two 
real-valued quantum fields (y9i(xi,t) and {p2{^i,t) as 



V5(Xi,t) = --^(V2l(xi,t) +2(y92(Xi,t)) 



^\^i,t) = -^(¥?l(Xi,t) -2(/?2(Xi,t)) 



one may separate the action as follows 

S = 3"^ + ^(1) + ^(2) + ^(3) ^ ^(4) 

u 



= N,J dt 



Jzq + fi — UuriQ 



j rft^(/?i(Xi,t) 



(16) 

(17) 
(18) 

(19) 



i a,b=l,2 



J r 

+ 27^*^ ^ </?a(Xi,t)¥?a(Xj,t) 



(20) 



<ij> a=l,2 
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99i(Xi,t)v9^(Xi,t) + V??(Xi,t) 



V9i(xi, t) + 9?2(xi, t) + 2v32(xi, t)ipl{:>Ci, t) 



(21) 
(22) 



In ( 120|) 5a;, is the antisymmetric tensor with £12 = 1, £21 = —1, and 

Xi = —fi + 3[/z/no 



(23) 



For a homogenous system the condensate is uniform and it is convenient to decompose 



the fluctuations into a Fourier series 

1 



22 



23| 



¥'a(Xj,t) 



E 



2i7rj 
AT., ' 



(24) 



where q = {qi, q2 ■ ■ ■ qd} with running from 1 to A^^ — 1 is an integer- valued vector field 
associated with all wave vectors in the BriouUin zone: q = 27rq/a, and 



jyl E = jyd E E • • • E • 



(25) 



The q = mode, i.e. the Goldstone mode, is omitted from the sum, to achieve orthogonality 
between the condensate and noncondensed modes. In momentum space the quadratic term 
S^'^^ as follows: 



q,q' V / 



(26) 



Mii(q,a;,q',a;') = -[Xi + e{ci) - J zq]5{uj + uj')5^^^^ , Mi2(q, q, w') = (27) 
M22(q,c<;,q,a;') = -[X2 + e{q) - J Zo]5{uj + uj')6^-c^ , M2i(q, w, q', w') = (28) 

with zq being the number of nearest neighbors. From this we extract the Fourier transfor- 
mation of the propagator of the fields y^i, and as the 2x2 matrix: 

/ 



G(a;,q) 



^2 _ + it 



X2 + £:(q) - Jzq -iu 
iuj Xi + e(q) — Jzq 



(29) 



where 



8{(i) = J(Xi + £(q) - Jzo){X2 + £(q) - Jzo 



£(q) = 2J[d-Y.cos{2T:qjN, 



(30) 



In coordinate space for a regular lattice the propagator is translational invariant 

Gab(xi,t;Xj,t') = Gabi^i - Xj,t - t') = (99a(xi,t)(^fc(xj,t')) 

Note that, in deriving (!26|) -( l30|) . the following relations have been used: 

"i27r,. 



H exp 

<mj> 



U • q - m • pj 



a=l 



I] exp 



q-p) 



III. THE EFFECTIVE POTENTIAL IN TWO-LOOP APPROXIMATION 



(31) 



(32) 



To organize the quantum corrections in a two-loop expansion, we separate the terms in 
the action f lT7|) into a free part and interaction parts following Jackiws pioneering work 20 1 



S 

Scl 



5i 



int 



'S'cl + 'S'fi-ee + 'S'int 

= Ns J dt |/iz/no + JzQvriQ — ij^^'^'^- 

^^V'a(Xi,t)M„6(Xi,t;Xj,t)v9f,(Xj,t) 
y (it^£int(v2l(Xi,t),V52(Xi,t)) 



(33) 
(34) 

(35) 
(36) 



(37) 



Ant(9'l(xi,t),¥J2(xi,t)) = V3[(/?l(xi,t)(/?|(xi,i) + 99f(xi,f)] 

+W4 \^\ (xi , t) + (/'I (xi , t) + 2v9f (xi , f ) (Xi , t)] = £3 + £4 

where 2x2 matrix Mah is given by Eqs. (1271) . (!28|) , fs = — t/ ^Jvnoj^, 

The perturbative framework is based on the propagator Gab{k,uj) given in ( l29l) . The 
effective potential V can be evaluated by the Eq. ( fT3l) . where the only connected, irreducible 
diagrams in the partition function Z = J VipiVip2GW {'^S{ipi,ip2)) should be taken into 
account. The grand thermodynamic potential i.e. free energy, Q{nQ,fi), corresponds to the 
minimum of V(no,/i), such that no is a solution of the equation dVijiQ, ji) /driQ = |24|. 



8 



Now using (1M|) - (136|) and making expansion by £int one can represent Z as follows: 

= e^^o / V^iVip2ei'P'^^''^''^^ 1 1 + ^ Ei / dtCi^tivi (xj, t), v^aCx;, t)) 

+ J dtdt'Cint {(Pl (Xi ,t),ip2 (Xi , t ) )£int (v^l (xj , t) , (/52 (xj , t ) ) | 

1 + i ^(AntCv^lCxi, t), (/'2(xi, t)))odi 



1.J 



(38) 



VDetG 



+ J (it(it'(Ant('/'l(xi,t),932(xi,t))£int(y'l(xj,f),(/?2(xj,t)))o > 

where we introduced the following notation 

^ ' G 

(i(v.„(xi, t), M^i,mo = A ( ^ f -r^) """^^ ' (39) 

V«c)ja(xi,t) ^c)J^,(xi,^)/ 

i=o 

suppressing the summation and integration signs over lattice sites and times t and t' in 
quandratic forms, for brevity. 

The classical contribution to V is given by factor exp (iSo) in f l38p 

Vo = ^ lne^^« = ^^[Uuno - 2(/. + Jzo)] (40) 

The one-loop contribution to the thermodynamic potential - Vn, can be obtained by using 
the free part of the action (13 5 p in f l38p . neglecting interaction terms: 

ViL = ^TVlnDetM = lnDetM(a;, q) (41) 



q 



where M(w,q) is given by (!28|) . One notices that the frequency sum, and with it the u 

n=oo 

integration, is divergent. In fact, to evaluate the frequency sum such as E ln(a^ + cu^), 

n=— oo 

with ujn = 2TTnT one differentiates it with respect to a and, after performing the summation 
over n, integrates it over a. This procedure gives an additional divergent constant which 
may be removed by an additive renormalization of the energy Therefore, in the case 
of optical lattices, where the momentum integration is performed within a finite volume 
there is no additional ultraviolet divergency coming from q integration, but there is an 
infinite constant coming from the frequency summation [26\. This divergent constant can 
be removed by subtraction from V the thermodynamic potential for the ideal gas 

ViT = Vii(f/) - v,l{u = 0) = ^ E^(q) - ^ E^(q) 

= ^E[^(q)-^(q) + /^ + ^^o], (42) 

^ q 



where we have used Eqs. (l23|) . (l30!l and performed integration by uj using formulas given in 
the Appendix. Further, for simphcity, we shall suppress the superscript in V{f^. 
The two-loop contributions to V are involved in second and third terms of as 

,•2 



V2L = ^ In < 1 + i ^ j {Cint)odt J dtdt' {CintjCint)o 



(43) 



The former includes £3(1^1,992) which does not contribute to Z, since it is in odd power of 
(fa, and hence: 



(>Cint)o = (A)o = V4{{(Pi)o + {^2)0 + '^{^1^2)0} 



(44) 



The same is true for {C2,{Va{'^i,t))'^4:{'^b{^i,t)))o coming from the third term of fl38|l . As to 
the term £4((y9a(xi, t))£4((y9a(xi, t)) it also should be omitted since its contribution is beyond 
two-loop corrections. Therefore 

V2L = ^Injl + z^(£4(<^i(xi,t),(^2(xi,t)))o 

+ ^ E / dtdt'iCsiM^i, t), ^2(Xi, t))£3(<^l(Xj, t), <^2(Xj, t)))o\ (45) 



The second term in the logarithm in Eq. fl45|) can be expressed in terms of propagator as 



(£4)0 = v^[3{Gl^{0) + Gl^m + 2Gii(0)G22(0) + 40^2 (0) 
where we used the following abbreviation x = (x, t) and the formulas 



((y9a(x)v9fe(x'))o 



Gab{x - X'), 

SGLIO), 

Gii(0)G22(0) + 26-^2(0), 



and introduced the notation 

Gab{0) = Gab{x,x) 



Note that 6*12(0) is the constant (see the Appendix) 



Gi2(0) 



E 

s q 



E[l] = -^21(0) 



(46) 



(47) 



(48) 



(49) 
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The third term, 

(£3^)0 = vl[{ipi{x)^l{x)^i{y)^l{y))^ + 2{ipi{x)^l{x)^^^^^^ (50) 
includes averages with six Lpa- These may be evaluated via Wick theorem to yield 

{ipl{x)ipi{y)ipl{y))o = 6Gn{x,y)Gl2{x,y) 

{ipi{x)ipl{x)ipi{y)ipl{y))o = 4:G22{x,y)Gi2{x,y)G2i{x,y) + 2Gl2{x,y)Gu{x,y). 

We have omitted one-particle reducible diagrams such as G22{0)Gii{x, y)Gu{0) . 
Now, using (gnD, (150])-([5I]) in (05]), we finally obtain: 

UNs 



V; 



2L 



3Gfi(0) + 3G^2(0) + 2G'n(0)G'22(0) + 4^^2(0) 



2T 



1.J 



J2 / c/trft'[G22(xi,t;Xj,t')G'ii(xi,t';Xj,t') 



+ 3G?i(xi, t; Xj, t') + 6G'n(xi, t; xj, tO^^lxi, t; xj, t' 
+ 2Gi2(xi,)f:;Xj,t')G'2i(xi,t;Xj,t')G'22(xi,t; ,Xj,t' 



(51) 



(52) 



The two-loop diagrams that contribute the thermodynamic potential are shown in Fig. 1. 





(a) 



m 





m 



(9) 



FIG. 1: Vacuum diagrams in a two-loop approximation. The solid and dashed lines correspond to 
Gil and G22 respectively, while the mixed line corresponds to G12 (or G2i)- 
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We now pass to momentum space, and perform integrations over energy variables u to 
obtain following analytic expression (see Appendix): 

Vii^ (no, /i) = |iV, (34 + 34 + 2/10/20 - 1) , (53) 
(no, ^^) = -^^{h + 3/2 - 6/3 + 2/4) . (54) 
where following integrals are introduced 

/io(no, ^) = = ^22(0), 

/2o(no, ^) = ^) = ^"(0)' 

(-/2 + 3Uuno + £(qi))(-/2 + SUunp + £(q2))(-/2 + ?7;vno + g(q3)) 
U^O,/^J - 2.q,^q, £:(qi)£:(q2)^:(q3)(^:(qi) +^:(q2) +^:(q3)) 

(55) 

(-/2 + ;7i/no + £(qi))(-/2 + Uunp + £(q2))(-/2 + [/t/no + £(q3)) 
2ino,/xj - £:(qi)£:(q2)£:(q3)(£:(qi) +£:(q2) +f (qs)) 

j( \ = V (-/2 + ?7i^no + £(q3)) 

i3lno, /ij £:(q3)(£:(qi) + £:(q2) + £:(q3)) ' 

X ^ {-fi- + 3Ui^no + ejqs)) 

i4(no,/.J - 2.,^^^^ £:(q3)(£:(qi) +£:(q2) +£:(q3)) 



and £:(q) = y (-yU + 3f/z/no + £:(q))Y'(-/i + f^^^^o + £:(q)), /i = /i - ^2:0, qs = qi - q2- 
Therefore the full effective potential in a two-loop approximation is given by 

V(/i,no) = Vo(/i,no) + ViL(/i,no) + V^^V^^o) + v£^(/i,no) (56) 

where Vo, Vil, , are given by equations (HOj) . (H2|) . (l53l) . (15^ respectively. Note 
that for the homogenous Bose gas. Eqs. fl53l) - fl56l) were calculated before by Braaten and 
Nieto H- 

IV. THE CONDENSATE FRACTION IN VPT 

To evaluate the condensate fraction no as an explicite function of U/J and u we shall use 



following strategy referred as a variational perturbation theory 



18l: 



1. With fixed values of input parameters introduce an auxiliary parameter, loop counter, 
r] {ri=l at the end of calculations) to represent V in Eq. fj56|) as: 

V(/i, no) = Vo(Ai, no) + ?7ViL(/i, no) + r]^V2L{fJ^, no) (57) 
12 



with V2L(/i, no) = V^iV? ^o) + v£V, rio) 

2. Impose the extremahzation condition: 

dV{fi,no) 



Oiv') (58) 



drio 

and solve this equation with respect to uq. Let the solution of the equation is nQ{fi). 
Clearly the latter can be also represented in powers of t]: 

no{ij) = noo(/i) + moiilA + ^^«02(^) (59) 

with 

noi(^)- nL(/^,-oo) 



?^02(/U) 



Vo(At,noo) 

(60) 

'^oi(/")K'(/^''^oo) + 2V22,(/^,noo) + 2noi(/i)Vii(/i,noo) 



2V^(/x,noo) 

where the prime denotes the derivative with respect to wq, e.g. V(^(/i, -n-oo) = 
[c^ViL(/i,no)/i9no]|no=noo and noo is the solution to the equation Vo(/i,no) = 0. 

3. Inserting n^^^) back to the effective potential fl57|) determines the free energy of the 
system fi(/i) = V(fiO)/^) 

4. Introducing a variational parameter M as 

H = M + rr] (61) 

with the abbreviation 

r = ^ (62) 

V 

and inserting f l6T|) into fi(/i) reexpand this Q{M,fi,r) in powers of ?7 at fixed r. 

5. Reinserting back r from ( l62l) optimize Q{M,^) with respect to the variational param- 
eter M. This will fix function of the optimal M = Mopt, with 

Mopt = Uu- Jzo (63) 

6. Finally, inserting this fi into fISIJl) one finds an explicit expression for no as no = 
no{U/J,iy). 
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Below we consider each step in detail. First, taking partial derivative with respect to no 
from Eq. fl57|) one presents fl58|) as 

(9V(no,/i) 9Vo(no,/i) , 9ViL(no,/i) ai^V^L (no, /i) ai^V^S (no, /i) 



c?no 
aVo(no,/i) 
dno 

dno 
;(i) 



c?no 



dno 



V- 



dno 



dno 



<9no 







-Ns [lyjl - Up' 



no 



f/zy (2/i — 3f/ z/no — 2£:(q) 



(64) 
(65) 
(66) 



^V^L (^o,Ai) f/^z^ (/" - ^(q)) [(/U -4f/z/no - £:(q))/io(no, /i) - (/i - £(q))/2o(no, /i)] 



E 



£3(q) 



(67) 



where following relations are used 
9£(q) Uv 



dno 
dho 
dno 

dho 



(2/i - Wuno - 2£(q) + 2Jzo) 



Up ■r-^ (/i — £^(q) + Jzo){n — SUuno — e{q) + Jzo 



2N. 



E 

s q 



£3(q) 



f/z/ (/i - £(q) + Jzo){^i - f/z/no - £(q) + J;zo) 



E 



dno 2N, ^ S^iq) 

In Eqs. ( 167|) 9V2/,V^^o has a long expression and will be given later. Solving Eq. (|6 
iteratively gives Eq. (15^ with 
/X + Jzo 

(^ + Jzo + 2e(q)) 



(68) 
(69) 
(70) 



nooif^) 
noiifJ') 



uU 
1 



2u 



no2iP') 



(3/2o(/i) + /io(Ai)) 

(2) 

1 91^2; ("-O'/") 



1 



E 



2^:^(q) 



5no 



1 

"0 = "-00 1 

^ 2[//2o(/x)g(q)(/x + J^o) ^ U{is + J zo)Hho{f^) - /20(/i)) 



C/e2(q)(/io(//)+/2o(/^)) 



(71) 



where 



/io(/i) = /io(no,/i)| 
^20 (yu) = ^2o(%,At)| 



"0=noo 



2iV, 



1 ^ 2/i + 2 j2;o + £(q) 



'^0=noo 



2iv,^^,,(q) 



(72) 



In this step the Goldstone boson dispersion is correctly achieved: 



SM) = V^(q) V^(q) + + 2 Jzo 



(73) 
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Now inserting (l59!l . (17T]) into (156|) one gets as a function of /i as 
n{fi) = V{fi, no) = Qoifi) + r]ni{fi) + rf'^2{p) 



Nsii^ + Jzof 
2U 



(74) 



02(//) = Vf2(/^,noo(/^)) + ^[2IM - 4/io(;U)/2o(^) - QI'M - 1]. 

Performed one more step of VPT we finally obtain as an explicit function of the parameters 
U, J, v. 



UN. 



Hq = Uu - Jzq, 



2N. 



fo(q) 



^^o(q) 



U ( /20B + 2 



+ 



UjhoB - hoB? 



, UHhoB+l20B-l)^e^{ci) 



and also the normal fraction, ni 



no as 



(75) 



E 



g(q) + [7;^ 
^^o(q) 



(76) 



n 



2L 



1 5V. 



(2) 



2L 



NMv"^ dm 



1 dV. 



(2) 
2L 



'^0='^00 



C/ [(/lOB - I20b){2U^i^^ - g'(q)) + ^t^£(q)(2/20i? - 1)] 



(77) 



E 



In Eqs. ([75]), ((TTj) £^o(q), Aos and /20B are given by 



hoB 



1 ^ 2Ui^ + g(q) 



(78) 



g(q) 

^^o(q)' 



Now we compare present approximation with Gutzwiller's. 
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n Gutzwiller approach the phonon dispersion for small q is quadraric in wave number 



12| rather than linear given in present approximation by Eq. ((73 



As it is seen from Eq.s f l76|) and (1771) in Bogoliubov type approximations the uncon- 
densed particles have momentum distribution Ug =< ajog > varying as for large 



momentum 



Il2|. 



29|, while in Gutzwiller approach this distribution is independent of q 



V. GROUND STATE ENERGY 



The ground state energy of the system at zero temperature can be determined as 

E = Q{fi) + fiN, 
where fl{fi) in Eq. (IHOj) can be rewritten as follows 

n{u, J, u) = no{u, J, v) + VLi{u, J, v) + n2{u, j, v), 



uNy 



n2{U,J,iy) = n^^l{U,J,iy) + 
where fl^l is given by 

^]g(f/,J,z/) = Vg(no = l,/i = /io) 



C/iV, (2/2q^ - 4/10B/20B - 6/|oB - 1) , Ns{fii-2Uufi2 



+ 



2U 



E 



4Ar2 

s qi,q2 



(79) 



(80) 



(81) 



£o(l)^o(2)^o(3)^oT ^o( 1)^0 (2)^0 (3)^oT 

with £o{q) given in (JTHl), and ei = e^^, £o{^) = ^^o(qi), £ot = ^^o(l) + £0(2) + £o{3). 

After some algebraic manipulations one obtains for the energy per particle E/N following 
expression 



N 8iy 2Uu 4z/ 

5:[£o(q)-^(q)]. 



2N,u q 



(82) 



Here the energy of an "ideal gas " (when [/ = in Bose- Hubbard Hamiltonian) has been 
subtracted . 
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VI. RESULTS AND DISCUSSIONS 



Firstly we discuss the condensate fraction, no vs U/J. In Fig. 2a it is presented in one- 
and two-loop approximations, (dashed and solid curves respectively) for the filling factor 
u = 1 and D = 3. It is seen that in the one loop approximation no can not reach zero within 
moderate values of U/J. More precisely no[one loop] = aX U/J = 81.2 . On the other 
hand, two-loop contributions coming from the diagrams in Fig. 1 are too large: quantum 
phase transition occurs at f//J ~ 6. Unfortunately this is rather far from the experimental 
value: no = at f//J ^ 29.34 as pointed out in the Introduction. It is seen from Fig. 2b that 
in Gutzwiller approach no reaches zero at [//J ~ 34.8 30|. Note that the similar behavior 
of no vs f// J with exactly the has been found by Stoof et al. in decouphng 



approximation in the second order perturbation theory 



16|. 





U/J 



20 
U/J 



a) 



b) 



FIG. 2: (Color online) The superfluid fraction no as a function of U/J for u = 1, D = 3. a)ln 
one (dashed line) and two loop approximations (solid line); b) Here the dashed line was obtained 
in Gutzwiller approache while the solid line in the present one. 

The superfluid fraction for two values oi u , u = 1 (dotted line) and u = 2 (solid line) 
is shown in Fig. 3a and Fig. 3b for D = 3 and D = 1 respectively. It is seen that the 
critical value of U/ J as well as a whole no(f/ / J, v) are not so sensitive to the filling factor. 
The fact that the superfluid fraction does not crucially depend on z/ has been observed 



also in Bogoliubov [9,] as well as HFB 



171] approximations. This is in contradiction with 
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2 
U/J 

a) 

FIG. 3: (Color online) The superfluid fraction as a function of U/J for u = 1 (dashed line) and 
V = 2> (solid line) for a) D = 3 and b) D = 1 in a two- loop approximation. 



10| where the dependence is rather 



the prediction by Gutzwiller single site approximation 
strong: 

Kcrit = ^o[\/^ + VlT^? = 2DyD+ VI + (83) 

Note that, although Eq. fl83l) gives a nice value for u = 1, Kent = 34.8, it can not be 
considered as an absolute truth since , besides it's drawbacks, outlined above, it takes into 
account the lattice dimensionality in a rather simple way. 

On the other hand as it is seen from Fig. 3b, for D = 1 the quantum phase transi- 
tion, which , more strictly speaking, is a Berezinskii- Kosterlitz- Thouless transition, occurs 
around U/J = 4. This is in good agreement with Monte- Carlo predictions 6[. Similar 
results for = 1 have been obtained by Danshita and Naidon in their time - evolving block 



decimation (TEBD) method 23|. However, note that TEBD method takes several days of 
computer calculations , while present approach does several minutes. In our calculations 
we used A^^ = 60, = uNg, that is we considered finite size systems. This explains the 
smoothness of no{U/J) in Figs. 3a and 3b. 

The ground state energy per particle E/N vs U/J in units Jzq in one (solid line) and 
two (dashed line ) loops is presented in Fig. 4. It is normalized such that the appropriate 
energy for the ideal case {U = in the Bose-Hubbard Hamiltonian) is set to zero. It is seen 
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4 

U/J 



8 



FIG. 4: (Color online) The energy per atom in units Jzq in one (solid line) and two-loop (dashed 
line) approximations for v = 1 for D = 2>. 

that quantum corrections due to diagrams in Fig.l are not significant for small U/J< 1. 
The dependence oi E/N on filling factor v is illustrated in Figs. 5a, 5b. It is seen that E is 
more sensitive to u than uq due to the leading term ( the first term in Eq. fl82p ) depending 
on V explicitly. 

VII. SUMMARY AND CONCLUSIONS 

We have developed a field theoretical approach in terms of path integral formalism to 
calculate the second-order quantum corrections to the energy density as well as to the 
superfiuid fraction in cubic optical lattices. Instead of using the standard formalism with 
complex field operatorsof condensed-matter literature, we find it more convenient to use two 
real fields. The thermodynamics of the system is deduced from the effective potential V, 
whose minimum gives free energy Vt. 

The superfiuid fraction, tiq , goes to zero at [// J ~ 6 for z/ = 1, 2, 3, and this is interpreted 
as a quantum phase transition from the superfiuid to the Mott insulator phase. For D = 1, 
we have found a good description of the transition. Unfortunately, for D = 2 and D = 3 the 
critical values for the parameters are rather far from the experiment: K^^f^{D = 2) = 16.8 
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4 8 4 

U/J U/J 



a) b) 
FIG. 5: (Color online) The energy per atom in units Jzq for various values of the filling parameter 
V for a) D = 3 and h) D = 1. 

and K'^f^^{D = 3) = 29.34, for u = 1. It appears that a more relaible value for Kent for 
D = 2,3 can only be reached by going beyond the present two-loop approximation. We 
expect that higher-order quantum corrections, for example post-Gaussian approximation 
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31| . will improve the situation, but they are hard to calculate. 
Thus we have shown that going beyond the Bogoliubov approximation employed by Stoof 
et al. 16|, one finds a quantum phase transition from a superfiuid to a Mott insulator state. 
Within a two-loop approximation we have derived explicit expression for the ground state 
energy of the optical lattice. 

Acknowledgment s 

We acknowledge support of the Volkswagen Foundation. AR is also indebted to the 
DAAD for partial support and to V. Yukalov and A. Pelster for useful discussions. 



20 



Appendix 

In present work all the calculations are carried out in real time. Loop integrals are taken 
over real energies u and over three dimensional quasimomentum k which pertains to the 
Brillouin zone —-nja < < / a. So, three or six dimensional integrals, presenting in one or 
two-loop calculations are finite and may be evaluated numerically by using Monte - Carlo 
methods. 

The integrals over uj are evaluated using contour integration. Some energy integrals 
needed for one- and two-loop calculations can be easily evaluated directly by using residue 
formulas: 



duj 



■oo 
-l-oo 

oo 
+00 

00 



27r (^2 _ 


- f 2 + ie) 


doj 


1 


2tt (a;2 - 




duj 




27r (cj2 _ 




1 





~2£ 

i 

i 



(84) 
(85) 
(86) 



X) 47r2 [uf - £f + ie][ujj - £l + ie][{ui + UJ2)^ - £l + ie] 4£:i£:2f3(fi + .fa + ^3 

dljJld0J2 UJ-\U}2 1 



cjiwa _ 
4^2 [^2 _ g2 + [^2 _ g2 ^ [(^^ +0J2y- £l + ie] ~ 4£^{£r+£^T£3) 

In the last two integrals £1 = £2 = £{^,2), and £^3 = £^(qi + qa)- 

The integral 

f duj 

^12 (q) = 



(87) 
(88) 



2^{u^ -£^{ci)+ie) ^^^^ 
needed for 6*12(0) = — (V^s) Z]q A2(q) should be considered more carefully. To evaluate it 
we use following formula given in the literature 26 1 



' n.= — r\ 



J7->0 



a(e/^» - 1) 



(90) 



where Un = 27mT, /3 = 1/T. The zero temperature limit, T — )■ 0, of fl90l) leads to 

^i2(q) = -I (91) 



so that Gi2(0) = {i/2Ns) J2q[^]- This constant enters into the evaluation of the constant 
m ~ (f^V), and produces a term —1 in the square brackets of Eq. (1761) . In a homogeneous 
Bose gas, such a constant term can be ignored. But here, on an optical lattice, it becomes 
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significant, so that in the evaluation of trace log term in Eq (HT]) , ^ it ^ must be taken into 
account properly. How to do that has been shown in the textbook [32]. Strictly speaking, 
the integral 

m = J^Hu'-£') (92) 

appearing in the trace log is divergent. To evaluate it, one may differentiate fl92|) with 
respect to S'^: 



(94) 



dS'^ J 27r(w2_^2) 

and use ( l8^ to obtain 

d£^ ~ 28 
Integrating this once £^ gives 

L{£) = [ ^ ln(a;^ - S"^) = i£ + constant. (95) 
J zvr 

Using the method of Ref. |32| we obtain the result of Section III where the constant leads 
to a term term — 1 in ni (see Eq. [76]) . 
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